home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / scheme / scm / jackal1a.lha / jacal / vect.scm < prev   
Encoding:
Text File  |  1992-12-24  |  9.8 KB  |  383 lines
  1. ;;; JACAL: Symbolic Mathematics System.        -*-scheme-*-
  2. ;;; Copyright 1989, 1990, 1991, 1992 Aubrey Jaffer.
  3. ;;; See the file "COPYING" for terms applying to this program.
  4.  
  5. (define (bunch_norm x)
  6.   (if (null? (cdr x))
  7.       (car x)
  8.     x))
  9.  
  10. (define (infinite-list-of . elts)
  11.   (let ((lst (copy-list elts)))
  12.     (nconc lst lst)))
  13.  
  14. (define (poly_elim poleqns vars)
  15.   (do ((vs vars (cdr vs)) (polys poleqns) (poly #f))
  16.       ((null? vs) polys)
  17.     (do ((var (car vs))
  18.      (pl polys (if (null? pl) 
  19.                (math-error "not enough equations" poleqns vars)
  20.                (cdr pl)))
  21.      (npl '() (cons (car pl) npl)))
  22.     ((poly_find-var? (car pl) var)
  23.      (set! poly (promote var (car pl)))
  24.      (do ((pls (cdr pl) (cdr pls)))
  25.          ((null? pls) (set! polys npl))
  26.        (if (bunch? (car pls)) (math-error "elim bunch?" (car pls)))
  27.        (set! npl (cons (poly_resultant poly (car pls) var)
  28.                npl))))
  29.       (if (bunch? (car pl)) (math-error "elim bunch?" (car pl))))))
  30.           
  31. ;;; This tries to solve the equations no matter what is involved.
  32. ;;; It will eliminate variables in vectors of equations.
  33. (define (eliminate eqns vars)
  34.   (bunch_norm
  35.    (if (some bunch? eqns)
  36.        (apply map
  37.           (lambda arglist (eliminate arglist vars))
  38.           (map (lambda (x)
  39.              (if (bunch? x) x (infinite-list-of x)))
  40.            eqns))
  41.      (poly_elim eqns vars))))
  42.  
  43. (define (copymatrix m)
  44.   (if (matrix? m)
  45.       (map copy-list m)
  46.       (math-error "Not a matrix: " m)))
  47.  
  48. (define (row? a)
  49.   (and (bunch? a)
  50.        (notevery bunch? a)))
  51.  
  52. (define (matrix? a)
  53.   (and (bunch? a)
  54.        (bunch? (car a))
  55.        (let ((len (length (car a))))
  56.      (and (every (lambda (r) (and (bunch? r) (= (length r) len)))
  57.              (cdr a))
  58.           len))))
  59.  
  60. (define (matrix . lst)
  61.   (if (matrix? lst) lst (matrix lst)))
  62.  
  63. (define (butnth n lst)
  64.   (cond ((null? lst) (math-error "Coordinate out of range: " n " " lst))
  65.     ((zero? n) (cdr lst))
  66.     (else (cons (car lst)
  67.             (butnth (+ -1 n) (cdr lst))))))
  68.  
  69. (define (mtrx_minor m i j)
  70.   (butnth (+ -1  i)
  71.       (map (lambda (x) (butnth (+ -1 j) x))
  72.            m)))
  73.  
  74. (define (transpose m)
  75.   (cond ((matrix? m)
  76.      (bunch_norm (apply map list m)))
  77.     ((bunch? m) (map list m))
  78.     (else m)))
  79.  
  80. (define (elim_test)
  81.   (define a (sexp->var 'A))
  82.   (define x (sexp->var 'X))
  83.   (define y (sexp->var 'Y))
  84.   (test (list (list a 0 0 124 81 11 3 45))
  85.     poly_elim
  86.     (list (list y (list x (list a 0 0 2) (list a 0 1)) 1)
  87.           (list y (list x (list a 5 1) 0 -1) 0 1)
  88.           (list y (list x (list a -1 3) 5) -1))
  89.     (list x y)))
  90.  
  91. (define (mtrx_square? m)
  92.   (if (not (matrix? m))
  93.       #f
  94.     (let ((rows (length m)))
  95.       (every (lambda (r) (= (length r) rows)) m))))
  96.  
  97. (define (mtrx_genmatrix fun i2 j2 i1 j1)
  98.   (do ((i i2 (+ -1 i))
  99.        (mat '()
  100.         (cons
  101.          (do ((j j2 (+ -1 j))
  102.           (row '() (cons (app* fun i j) row)))
  103.          ((< j j1) row))
  104.          mat)))
  105.       ((< i i1)
  106.        (bunch_norm mat))))
  107.  
  108. ;;;Create a matrix of indeterminates
  109.  
  110. (define (nmatrix var . rst)
  111.   (apply genmatrix
  112.      (lambda subs (apply make-subscripted-var var subs))
  113.      rst))
  114.  
  115. (define (dotproduct a b)
  116.   (let ((cols (matrix? a))
  117.     (colbs (matrix? b)))
  118.     (cond ((or (and cols (bunch? b)) (and colbs (bunch? a)))
  119.        (if (not (and (eqv? cols colbs) (eqv? (length a) (length b))))
  120.            (math-error "dotproduct of unequal size matrices: " a
  121.                " " b))
  122.        (reduce (lambda (x y) (app* _@1+@2 x y))
  123.            (reduce (lambda (x y) (app* _@1+@2 x y))
  124.                (map (lambda (x y) (app* _@1*@2 x y))
  125.                 a b))))
  126.       ((and (row? a) (row? b))
  127.        (reduce (lambda (x y) (app* _@1+@2 x y))
  128.            (map (lambda (x y) (app* _@1*@2 x y))
  129.             a b)))
  130.       ((bunch? a) (map (lambda (x) (app* _@1*@2 x b))
  131.                a))
  132.       ((bunch? b) (map (lambda (x) (app* _@1*@2 a x))
  133.                b))
  134.       (else (app* _@1*@2 a b)))))
  135.  
  136. (define (ncmult a b)
  137.   (let ((cols (matrix? a))
  138.     (colbs (matrix? b)))
  139.     (cond ((or (and cols (bunch? b)) (and colbs (bunch? a)))
  140.        (if (not colbs) (set! b (list b)))
  141.        (if (not (= (or cols (length a)) (length b)))
  142.            (math-error "matrix product of unequal size matrices: " a
  143.                " " b))
  144.        (or cols (set! a (matrix a)))
  145.        (bunch_norm
  146.         (map (lambda (arow)
  147.            (apply map
  148.               (lambda bcol
  149.                 (dotproduct bcol arow))
  150.               b))
  151.          a)))
  152.       (else (deferop 'ncmult a b)))))
  153.  
  154. (define (mtrx_expt a pow)
  155.   (cond ((zero? pow) 1)
  156.     ((one? pow) a)
  157.     ((negative? pow) (mtrx_expt (mtrx_inverse a) (- pow)))
  158.            (else
  159.      (ipow-by-squaring a pow 1 dotproduct))))
  160.  
  161. (define (crossproduct a b)
  162.   (if (and (row? a) (row? b))
  163.       (cond ((not (= (length a) (length b)))
  164.          (math-error "Cross product of unequal length vectors: " a
  165.              " " b))
  166.         ((= 2 (length a))
  167.          (app* _@1-@2
  168.            (app* _@1*@2 (car a) (cadr b))
  169.            (app* _@1*@2 (car b) (cadr a))))
  170.         ((= 3 (length a))
  171.          (list
  172.           (app* _@1-@2
  173.             (app* _@1*@2 (cadr a) (caddr b))
  174.             (app* _@1*@2 (cadr b) (caddr a)))
  175.           (app* _@1-@2
  176.             (app* _@1*@2 (caddr a) (car b))
  177.             (app* _@1*@2 (caddr b) (car a)))
  178.           (app* _@1-@2
  179.             (app* _@1*@2 (car a) (cadr b))
  180.             (app* _@1*@2 (car b) (cadr a)))))
  181.         (else (math-error "Cross product of funny length vectors: " a
  182.                   " " b)))
  183.       (dotproduct a b)))
  184.  
  185. (define (mtrx_scalarmatrix n x)
  186.   (do ((i 1 (+ 1 i))
  187.        (rows (list (nconc (make-list (+ -1 n) 0) (list x)))
  188.          (cons (append (cdar rows) (list 0)) rows)))
  189.       ((>= i n) rows)))
  190.  
  191. (define (mtrx_diagmatrix diag)
  192.   (set! diag (reverse diag))
  193.   (do ((i (+ -1 (length diag)) (+ -1 i))
  194.        (j 0 (+ 1 j))
  195.        (diag diag (cdr diag))
  196.        (rows '()
  197.          (cons (nconc (make-list i 0) (list (car diag)) (make-list j 0))
  198.            rows)))
  199.       ((null? diag) rows)))
  200.  
  201. ;;; Something is subtly worng here.
  202. (define (determinant m)
  203.   (if (not (mtrx_square? m)) (eval-error "determinant of non-square matrix" m))
  204.   (letrec
  205.       ((x11 1)
  206.        (ds (lambda (m)
  207.          (cond ((null? m) 0)
  208.            ((expr_0? (caar m))
  209.             (set! x11 (- x11))
  210.             (cons (cdar m) (ds (cdr m))))
  211.            (else
  212.             (set! x11
  213.               (if (negative? x11) (app* _-@1 (caar m)) (caar m)))
  214.             (map (lambda (ro)
  215.                (if (expr_0? (car ro))
  216.                    (cdr ro)
  217.                    (app* _@1-@2*@3
  218.                      (cdr ro)
  219.                      (cdar m)
  220.                      (app* _@1/@2 (car ro) (caar m)))))
  221.              (cdr m)))))))
  222.     (cond ((null? (cdr m)) (caar m))
  223.       ((null? (cddr m)) (crossproduct (car m) (cadr m)))
  224.       (else (let ((subdet (determinant (ds m))))
  225.           (app* _@1*@2 x11 subdet))))))
  226.  
  227. (define (charpoly m var)
  228.   (determinant (app* _@1-@2
  229.              m
  230.              (mtrx_scalarmatrix (length m) var))))
  231.  
  232. (define (cofactor
  233.      m i j)
  234.   (let ((det (determinant (mtrx_minor m i j))))
  235.     (if (odd? (+ i j))
  236.     (app* _-@1 det)
  237.     det)))
  238.  
  239. (define (coefmatrix eqns vars)
  240.   (map
  241.    (lambda (eq)
  242.      (map (lambda (var)
  243.         (let ((c (poly_coeff eq var 1)))
  244.           (set! eq (poly_coeff eq var 0))
  245.           c))
  246.       vars))
  247.    eqns))
  248.  
  249. (define (augcoefmatrix eqns vars)
  250.   (map
  251.    (lambda (eq)
  252.      (nconc
  253.       (map (lambda (var)
  254.          (let ((value (poly_coeff eq var 1)))
  255.            (set! eq (poly_coeff eq var 0))
  256.            value))
  257.        vars)
  258.       (list eq)))
  259.    eqns))
  260.  
  261. (define (echelon m) m)
  262.  
  263. ;;; This is not correct
  264. (define (triangularize m)
  265.   (let* ((cols (matrix? m))
  266.      (n (length m))
  267.      (c (make-list n -1)))
  268.     (cond ((null? cols) (set! m (list m))
  269.        (set! n 1))
  270.       ((> cols n) (set! m (copymatrix m)))
  271.       (else (set! m (transpose m))
  272.         (set! n cols)))
  273.     (do ((k 0 (+ 1 k)))
  274.     ((>= k n) m)
  275.     (let* ((j
  276.         (do ((j 0 (+ 1 j)))
  277.             ((>= j n) #f)
  278.             (if (and (< (list-ref c j) 0)
  279.                  (not (poly_0? (list-ref (list-ref m k) j))))
  280.             (return j))))
  281.            (rowj (and j (list-ref m j))))
  282.       (cond (j
  283.          (set! rowj (app* _-@1/@2
  284.                   rowj
  285.                   (list-ref rowj k)))
  286.          (set-car! (list-tail m j) rowj)
  287.          (do ((restrows m (cdr restrows)))
  288.              ((null? restrows))
  289.              (if (not (eq? (car restrows) rowj))
  290.              (set-car! restrows
  291.                    (app* _@1*@2+@3
  292.                      rowj
  293.                      (list-ref (car restrows) k)
  294.                      (car restrows)))))
  295.          (set-car! (list-tail c j) k))
  296.         (else #f))))))
  297.  
  298. ;;; This algorithm taken from:
  299. ;;; Knuth, D. E.,
  300. ;;; The Art Of Computer Programming, Vol. 2: Seminumerical Algorithms,
  301. ;;; Addison Wesley, Reading, MA 1969.
  302. ;;; pp 425-426
  303. (define (rank m)
  304.   (let* ((cols (matrix? m))
  305.      (n (length m))
  306.      (c (make-list n -1))
  307.      (r 0))
  308.     (cond ((null? cols) (set! m (list m))
  309.        (set! n 1))
  310.       ((> cols n) (set! m (copymatrix m)))
  311.       (else (set! m (transpose m))
  312.         (set! n cols)))
  313.     (do ((k 0 (+ 1 k)))
  314.     ((>= k n) (- n r))
  315.     (let* ((j
  316.         (do ((j 0 (+ 1 j)))
  317.             ((>= j n) #f)
  318.             (if (and (< (list-ref c j) 0)
  319.                  (not (poly_0? (list-ref (list-ref m k) j))))
  320.             (return j))))
  321.            (rowj (and j (list-ref m j))))
  322.       (cond (j
  323.          (set! rowj (app* _-@1/@2
  324.                   rowj
  325.                   (list-ref rowj k)))
  326.          (set-car! (list-tail m j) rowj)
  327.          (do ((restrows m (cdr restrows)))
  328.              ((null? restrows))
  329.              (if (not (eq? (car restrows) rowj))
  330.              (set-car! restrows
  331.                    (app* _@1*@2+@3
  332.                      rowj
  333.                      (list-ref (car restrows) k)
  334.                      (car restrows)))))
  335.          (set-car! (list-tail c j) k))
  336.         (else (set! r (+ 1 r))))))))
  337.  
  338. (define (matrix_test)
  339.   (let ((d2 (list (list 2
  340.             (list a 1 1)
  341.             (list b 0 -5))
  342.           (list (list a 0 1)
  343.             (list b 0 1)
  344.             (list c 0 1)))))
  345.     (test (list (list 1 (make-rat 1 2) (make-rat 1 3))
  346.         (list (make-rat 1 2) (make-rat 1 3) (make-rat 1 4))
  347.         (list (make-rat 1 3) (make-rat 1 4) (make-rat 1 5)))
  348.       mtrx_genmatrix
  349.       (make-rat 1 (list _@2 (list _@1 -1 1) -1)) 3 3 1 1)
  350.     (test d2
  351.       augcoefmatrix
  352.       (list (list y (list x (list b 0 5) -2) (list a -1 1))
  353.         (list y (list x (list c 0 1) (list a 0 1)) (list b 0 1)))
  354.       (list x y))
  355.     (test (list (list 1
  356.               (make-rat (list a 1 1) 2)
  357.               (make-rat (list b 0 5) 2))
  358.         (list 0
  359.               1
  360.               (make-rat (list c (list b 0 (list a 0 5)) -2)
  361.                 (list b (list a 0 -1 -1) 2))))
  362.       echelon
  363.       d2)
  364.     (test (list (list 2
  365.               (list a 1 1)
  366.               (list b 0 -5))
  367.         (list 0
  368.               (list b (list a 0 -1 1) 2)
  369.               (list c (list b 0 (list a 0 5)) 2)))
  370.       triangularize
  371.       d2)
  372.     (test 2
  373.       rank
  374.       d2)
  375.     (test (list y 10 -7 1)
  376.       charpoly
  377.       (list (list 3 1)
  378.         (list 2 4))
  379.       (list y))))
  380.  
  381. ;;;    Copyright 1989, 1990, 1991, 1992 Aubrey Jaffer.
  382.  
  383.